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We perform off-lattice, canonical ensemble molecular dynamics simulations of the self-assembly 
of long segmented copolymers consisting of alternating, tunably attractive and hydrophobic hinder 
domains, connected by hydrophilic linker chains whose length may be separately controlled. In such 
systems, the molecular design of the molecule directly determines the balance between energetic 
and entropic tendencies. We determine the structural phase diagram of this system, which shows 
collapsed states (dominated by the attractive linkers’ energies), swollen states (dominated by the 
random coil linkers’ entropies) as well as intermediate network hydrogel phases, where the long 
molecules exhibit partial collapse to a single moleeule network state. We present an analysis of the 
connectivity and spatial structure of this network phase, and relate its basic topology to mechanical 
properties, using a modified rubber elasticity model. We find that it is possible to optimize the 
mechanical performance by an appropriate choice of molecular design, which may point the way to 
novel synthetics that make optimal mechanical use of constituent polymers. 

PACS numbers: 82.20.Wt, 82.35.Jk, 81.16.Fg 


I. INTRODUCTION 

Biological materials consist mostly of networks of fil¬ 
amentous fibers. A prime example is the cytoskeleton, 
which - among its many other functions - governs me¬ 
chanical stability, the response to stress, and the motil¬ 
ity of the cell. Its mechanical properties, particularly in 
the non-linear regime are uniquely linked to its molecular 
architecture: an open meshworks of interconnected semi- 
flexible fibers mm- Mimicking these mechanical prop¬ 
erties in a synthetic material would offer important new 
possibilities both for cell research and health technologies 
such as drug delivery or cell growth, as well as for novel 
bio-insipred and biobased performance materials. In this 
paper, we investigate a self-assembly route to recreating 
the topology of such materials in synthetic associative 
polymers, and ask the question how these should be de¬ 
signed in order to spontaneously produce networks that 
present optimal mechanical performance. To this end 
we study, numerically, a self-assembling system of block 
copolymers that forms an effective cross-linked network. 

The self assembling system we focus on here is a multi¬ 
block amphiphile, consisting of many regular repeats of 
alternating hydrophilic and hydrophobic blocks. Previ¬ 
ous work, which we will briefly review in a moment, has 
identified such molecules as quite promising candidates 
for precise regulation of self-organized architecture, but 
so far experimental realizations have largely been lack¬ 
ing. Polysaccharides with hydrophobic substitutions - 
methylcellulose, in particular - possess the required al- 
ternatingly hydrophilic/hydrophobic backbone structure, 
but in more or less randomly distributed blocks [3]. Re¬ 
cently, however, a novel molecule was synthesized and 
studied [4] that does allow precise control over the spa¬ 
tial arrangement of the alternating blocks, and in fact 
also over the attractiveness of the individual hydropho- 
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FIG. 1: We simulate long copolymers consisting of alternating 
binder (yellow) and linker (blue) blocks. The binder blocks 
consist of repeating PEG units, the linker blocks are bisurea 
motifs. In synthesis, the length of both these blocks is sep¬ 
arately tuneable. In our simulations, we coarse grain these 
molecules into beads representing many atoms at once. The 
molecular design is fixed by specifying values for n, m and k. 


bic blocks. Its basic building blocks (shown in Fig. are 
macromolecules consisting of alternating binder (attrac¬ 
tive) and linker blocks, which are, respectively, bis-urea 
hydrophobic motifs (essentially identical to those dis¬ 
cussed in [5]), and polyethylene glycol (PEG) hydrophilic 
chains [4]. The objective of this work is to provide an ini¬ 
tial survey, based upon Molecular Dynamics simulation 
and elementary rubber elastic theory, of the expected mi¬ 
crostructure of the self-assembled states of this molecule, 
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as well as the mechanical properties of its gel/network 
phase. 

In our simulations, we coarse grain many repeat units 
into beads that interact via standard hard-core repulsive 
and - for the binder blocks - attractive potentials. 

In aqeous solvent the amphiphilic design of these 
molecules leads to aggregation and causes the long 
molecules to partially collapse into clusters that are con¬ 
nected via one or more links. Clearly, if the attraction 
is too weak the chains will adopt random coil configura¬ 
tions, while if the attraction is too large the molecules will 
collapse on themselves forming compact globular aggre¬ 
gates. In between, as we will show, there is an interme¬ 
diate regime where the system spontaneously aggregates 
into a network phase of connected clusters linked to each 
other. Here we present simulations in which we control 
the linker length and the interaction potential between 
the binders to tune the relative importance of energy 
and entropy in this self-assembly process, and charac¬ 
terize the resultant network in terms of cluster size and 
connectivity. We summarize these findings in a phase 
diagram, identifying the experimentally most promising 
regime for recovering biopolymer network-like architec¬ 
tures in self-assembling copolymers. We then analyze 
the mechanical properties of these networks, first in the 
context of a rubber-elastic theory for networks of vari¬ 
able functionality, then in direct numerical simulations. 
We are not the first to consider models of self assembling 
multiblock copolymers: Rooted in extensive work on the 
assembly and rheology of amphiphilic block copolymers 
(see, e.g., i and references therein) more recent work has 
also explicitly considered the multiblock design. Early 
scaling theory by Halperin [7] anticipated the emergence 
of flowerlike micellar structures, either single or multiple 
such structures connected to each other in the manner of 
pearls on a string, depending on the quality of the sol¬ 
vent. Subsequent efforts refined the conditions for 

formation and stability of the flower like micelles. Later 
numerical work unmi] confirmed that indeed such struc¬ 
tures may form in dilute systems. In [l2j[T3], the Monte 
Carlo (MC) technique on a lattice model was applied on 
the system to establish an equilibrium phase diagram, 
which features characteristic phases of isloated micelles, 
connected micelles, but also laminar and tubular phases. 
Similar results were reported in nnun] in yet more MC 
simulations. More recently, the first dynamical simula¬ 
tions (Brownian Dynamics, BD) were presented, varying 
systematically the solvent quality. Again, the familiar 
structural phases were reported m- The picture that 
emerges is a very rich self-assembling system, which - pro¬ 
vided one has some synthetic control over the moelcular 
architecture - provides a direct control over mesoscopic 
organization of a hydrogel. A question that has received 
very little attention, so far, has been the implications of 
such structure for mechanical performance. Our ultimate 
ambition is to simulate dynamical rheology for these net¬ 
works, in conjunction with the self-assembly. We focus on 
those structural aspects of the hydrogel phase that govern 


mechanical response. This response is analyzed first in 
the context of simple rubber elastic models, and then di¬ 
rectly, using a computational implementation of the oscil¬ 
latory rheology protocol on the networks presented here. 
While - aside from the polysaccharides mentioned above 
- there are few experimental realizations to compare to, 
the molecular design we coarse grain here was, in fact, 
studied in previous experimental work [4] , which revealed 
a tendency to form strong hydrogels, whose rheological 
properties make them uniquely suited as injectable sub¬ 
strates for drug delivery. However, despite the detailed 
AFM and SAXS analysis of the self-assembled hydrogels, 
the question of the precise molecular structure of the gels 
remained unanswered. In this paper, we use molecular 
dynamical simulations to address this question, which 
enables us to determine a relation between molecular de¬ 
sign, gel structure, and mechanical properties that breaks 
ground for further rational design in these materials. 

This paper is organized as follows: Section [HI presents 
the experimental model system, the modeling environ¬ 
ment and approach we have adopted, and the relevant 
interaction potential we employ. In Section HI the sim¬ 
ulation protocol in coarse-graining the model system is 
presented. In Section [TV| we present data on the phase 
behavior resulting from the self-assembly process. Sec¬ 
tion |V| analyzes the topology of the network phase and 
presents results on the size and distribution of the con¬ 
nectivity of clusters. In the remainder of the paper we 
investigate the mechanical response of these networks, 
in Section \VJ\ we recall and apply a modified rubber 
elastic model, in Section [VH] we measure the rheological 
response of such networks numerically, using a computa¬ 
tional bulk rheology approach. This section is followed by 
qualitative result section. We then relate these findings 
to experimental observations. We present our conclusions 
for the relation between molecular design, supramolecu- 
lar network structure and the mechanical properties in 
our system. 


II. MODEL SYSTEM AND COARSE GRAINED 
MD SETUP 

As indicated in Fig. the multiblock copolymer 
we study here is a polymer consisting of repeat¬ 
ing diblock units. Its structure may be denoted as 
[[PEG]n[bisurea^]]/c; n, m and k are the repeat numbers 
of the PEG molecule, the bisurea motif, and the resul¬ 
tant diblock unit respectively. In typical experimental 
settings, n ^ 10^,m ^ 5,/c ^ 10 |4]. The PEG chain 
is hydrophilic and highly flexible which results in a 
tendency to swell in water due to entropy maximization. 
The bis-urea block, in contrast, is stiff (nonswellable) 
and hydrophobic, resulting in a generic tendency to 
aggregate in water. This tendency is further enhanced 
in by the ability of the ureas to hydrogen bond, leading 
to an increased stability of their aggregates. Thus, the 
ultimate behavior of this long copolymer is determined 
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by conflicting tendencies imparted by binders and 
linkers; depending on the dominating feature the system 
will lean towards collapsed (binder dominated) states 
or swollen (linker dominated) states. Interestingly, the 
geometry of the molecule may be used to directy affect 
this balance between energy and entropy. The longer the 
linker chain (n), the larger the entropic contribution; the 
longer the bisurea block (m) the stronger the hydropho- 
bicity/hydrogen bond-driven tendency to collapse. 



FIG. 2: Graph of the Lennard-Jones (LJ) potentials between 
binder beads (left). Linker beads, interacting with either 
other linker beads or binder beads do so with a LJ poten¬ 
tial truncated at r = a to approximate hard-core repulsion 
between beads(right). 

In our simulations, we are initially interested in the 
generic features of this class of systems. We therefore 
model a copolymer chain as a sequence of permanently 
attached beads of fixed (and equal) radius a. There are 
two distinct types of beads, and they are distinguished 
by their mutual interactions (see Fig. |^. All beads 
repel strongly at distances shorter than cr, to enforce 
hard-core repulsion and non-crossing of the copolymer 
(we address this in more detail below). Binders mutu¬ 
ally attract at distances larger than a, whereas linkers 
experience no interactions, either with other linkers or 
binders, beyond their hard-core radius. As we are ini¬ 
tially interested in open, meshworked phases we shall 
use as our two main control parameters in this system 
the number of binder beads the number of linker 
beads Ni, and the strength of the LJ attraction as mea¬ 
sured by the LJ well depth 5 . As the experimental 
molecules have either m = 4 or m = 6, the physical 
size of the binder region does not vary much, which is 
why we choose to fix = 3, reducing the number of 
free parameters. Simulated coarse grained molecules are 
denoted as B(A^ 5 )L(A^/), so that, for instance 531/24 de¬ 
notes a molecule whose repeating diblock motif consists 
of 3 binder beads and 24 linker beads. For all simula¬ 
tions reported here, we use 500 repeats of this diblock, 
so that each molecule in each simulation contains 1500 
binder beads. 

As was previously observed in lattice MC and BD sim¬ 
ulations [mini ns, the equilibrium phase diagram of 
such systems is already very rich. At low values of 5 , 
or very large linker regions, the entropy dominates and 



FIG. 3: Evolution of the self-assembly of a network during a 
typical MD run in the network phase. Starting from a random 
arrangement, the macromolecule gradually collapses on itself 
to form clusters of a well defined size, connected by one or 
more linker chains. Note, that these images are still of single 
(albeit very large) molecules. 


swollen phases are expected. For very high values of the 
LJ attraction, or for short linkers, energy dominates and 
complete collapse (similar to a polymer in poor solvent) 
is observed. In the intermediate regime, chains of mi¬ 
celles are reported, as are nonspherical aggregates. The 
objective of the present paper is to establish whether 
these structures are also present in off-lattice MD sim¬ 
ulations. As will show, most do indeed feature in our 
simulations, although a large region of the dynamical 
phase space is occupied by meshlike structures which are 
best termed flower like micellar networks, which appear 
to be dynamically arrested intermediates that are unable 
to fully equilibrate to a single collapsed state, see Fig. 

As these flower like micellar networks most closely re¬ 
semble the biopolymer gels whose topology we aim to 
recreate, we characterize the connectivity structure and 
relate it to mechanical response using modified rubber 
elastic theory. 

Our simulations are carried out using the LAMMPS 
molecular dynamics package m- We represent the 
long copolymers using a bead-spring model, under fixed 
boundary conditions, with two distinct types of beads to 
represent binder and linker segments. Adjacent beads 
interact via a harmonic bond potential of the form 
Ubondi'f^) = kb{r — lb)‘^, and a harmonic bending potential 
Uangie{^) = ka{0 — tt)^, is applied to each set of three 
neighboring binder beads (where r is the distance be¬ 
tween the centers of mass of pairs of beads, = 1.2 a is 
the equilibrium bond length, a the size (diameter) of the 
beads and 0 the angle formed by the two bonds that con¬ 
nect the middle bead to its two adjacent beads), and 
ka are the stretching and bending stiffnesses, which are 
both fixed at fairly high values: k^ = 200 kBTjcP' and 
ka = 400, to enforce inextensibility of the back¬ 
bone, and inflexibility of the binder region during the 
simulations. The bending stiffness is set to zero for linker 
beads, so that the linkers are represented as flexible poly¬ 
mers of molecular weight proportional to Ni. All binder 
beads that are not first or second nearest neighbors (i.e., 
those that are part of different binder domains) inter¬ 
act through a full Lennard-Jones (LJ) potential, whose 
attractive strength measured by 5 : 
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The interactions involving linker beads (i.e, linker- 
linker and binder-linker) are purely repulsive, and are 
modeled with a truncated LJ potential, (see Fig. 


U{r) 



if r < cr 
if r > cr. 


( 2 ) 


In the latter equation, e^ep measures the strength of 
the repulsion, which in our simulations is set equal to 
0.01 /cbT; r is still the center-to-center distance between 
the beads. The repulsive potential is used to ensure that 
the polymer obeys self-avoidance: We geometrically en¬ 
sure that the chain cannot cross itself by setting the bead- 
bead distance to 1.2cr, and we truncate the LJ potential 
at r = cr. The repulsive LJ potential takes care of the 
self avoidance of all beads, and to prohibit self-crossing 
of the polymer we choose an equilibrium bond length la 
of 1.2cr. This leaves a gap with an equilibrium size of 
0.2cr between neighbouring beads, and as one may see in 
Fig-i the bond length never fluctuates to lengths suffi¬ 
cient to leave a space through which another segment of 
polymer may pass. 


0.04 


0.03 

mio-) 

0.02 


0.01 


0.00 

1.15 1.20 1.25 


bond < CT V 2 



4/0- 


FIG. 4: Histogram of bond length distribution during one run. 
The spacing between two beads is never sufficiently large to 
allow other beads through, and hence the chain cannot self 
cross during the dynamics. 


may overlap slightly. As one may see from Fig. the 
exponent u rises quickly from 0.5 (the ideal chain value) 
to its self-avoiding value upon 0.6 for increasing values of 
^rep- We fix 5rep at 0.01 becausc at this value, we first see 
(within our numerical accuracy) the correct exponent of 
z/ = 0.6. 


III. SIMULATION PROTOCOL 

In our LAMMPS simulations, each simulation is re¬ 
peated from random initial conditions for 5 different sys¬ 
tems, which are then averaged over. Our protocol is 
as follows: We begin each simulation by an equilibra¬ 
tion process starting from a random initial geometry. 
To capture the effects of different attractive strengths 
(representing the different cohesive energies, which rises 
with m in the molecule) and different linker lengths, 
we simulate twelve types of coarse grained molecules: 
B3L3, B3L4, B3L5, B3L6, B3L9, B3L12, B3L18, B3L24, 
B3L45, B3L75, B3L120, B3L150, and each of these is 
simulated for twelve values of 5. After an initial random¬ 
ization, we turn on the attractive interactions and watch 
the system evolve. 



FIG. 5: Rg scaling for contour lengths of 20, 50,100, 200 and 
300 beads. The attractive potentials are switched off. The 
Flory exponent for beads of diameters of d = 0.5,1 and 2 is 
in agreement with theoretical prediction. 


IV. PHASE DIAGRAM 


To verify that indeed, this reproduces correct polymer 
scaling, we turn off the attractive interactions (setting 
5rep, as well as the bond-bending term) and record the 
radius of gyration as a function of the length of the poly¬ 
mer. As one may see in Fig. we recover the correct 
self-avoiding scaling Rg ^ with a Flory exponent u of 
about 0.6. This scaling is unaffected by the actual size 
of beads. The repulsive LJ potential [^exists between all 
beads, its value affects the effective radius of a bead. If 
we choose it to be too small, the beads become ’soft’ and 


What happens next obviously depends on the strength 
of the attractive interactions, and the geometry of the 
molecule. For low values of 5, much like swollen polymers 
in a good solvent, we see random coil configurations, 
determined by the entropy of the linkers. Upon in¬ 
creasing 5, the chains start to self-assemble and form 
clusters. We see the size of clusters grow with increasing 
5 but their growth rate is limited by the entropy of the 
linkers. For short linkers, collapse to a single cluster 
at high 5 values is observed. For even longer linkers. 
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FIG. 6: Scaling exponent v for different Srep- We have chosen 
£rep — 0.01 in all our simulations which as this is the value 
where correct real chain scaling is first observed. 

however, the average cluster size remains limited since 
the entropic penalty for packing the linkers into the 
corona of the aggregate rises strongly with their length. 
The intermediate regime is where we see a network of 
flowerlike micellar cores (small aggregates of binders), 
linked to other such clusters by one or more linker 
chains - see Fig. This phase diagram is in good 

agreement with a similar result reported in earlier lattice 
Monte Carlo simulations, though we obviously cannot 
reproduce the lattice-induced layered aggregate phase 
reported there. 

A second, more striking difference is that the network 
phase we see appears to be more prevalent in our MD 
simulations than in the earlier MC work [HIISIIIS]. We 
attribute this to kinetic trapping of the structures: while 
truly long relaxation might recover collapsed phases as 
true thermal equilibria, on the timescales that we are 
able to access in our MD sims the system tends to arrest 
in long lived, but possibly metastable, stationary states. 
In the following, we characterize the topology in this net¬ 
work phase. 


V. NETWORK STRUCTURE ANALYSIS 

Clearly, the region in the phase diagram that most 
closely resembles that of a crosslinked biopolymer mesh 
is the network phase. As is known from extensive pre¬ 
vious work [DlllIIIHlI], the topology of such a network 
(connectivity, crosslinker density) greatly influence the 
mechanical response. We will quantify the connectivity 
and cluster size of the networks that form in the interme¬ 
diate regime, as these parameters are expected to directly 
translate into effective parameters for the crosslinked net¬ 
work: We identify each cluster with a network node, 
whose cohesive energy is determined by the number of 



FIG. 7: The phase diagram of the self-assembled states of 
the macromolecule is divided into three distinct phases: large 
single clusters (violet), single molecule networks (green) and 
random coils (purple). The network phase is also loosely di¬ 
vided into three sub-phases: One with a broad distribution of 
cluster sizes (0), ane with much tighter distribution of cluster 
sizes (•), and one with a combination of clusters and isolated 
single beads (A). 


binders it contains, and we count the number of connect¬ 
ing links between such network nodes to determine the 
distribution of functionalities in the network. 



Ns 


FIG. 8: Frequency distribution of the step length Ns for B3L7, 
at e: = 1.5. The rapid decay justifies the approximation we 
will make, which is that all linkers are one step linkers. 

Our method to identify the clusters is as follows: Af¬ 
ter equilibration, we export the MD trajectories of all 
beads. To identify individual clusters, we sort all binders 
that are within dmax from each other in a group, dmax is 
the mean value of the longest distance between two cen¬ 
ter of mass (3.4cr) (the heart-to-heart distance between 
the central beads of two aligned binder blocks) and the 
shortest such distance (a) (i.e., the heart-to-heart dis¬ 
tance between the central beads of two parallel binder 
blocks). The average cluster size is then an average over 
the distribution of number of binders in each group. This 
computation is averaged over five separate runs. We also 
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calculate the functionality f of each cluster, defined as 
the number of linker chains that bridge to other clusters 
emanating from one particular cluster. In our simulations 
we count all of the chains that connect two separate clus¬ 
ters together. However, the length of these bridges can 
be one or more times the linker length, and the length 
of the linker determines its effective elasticity. To keep 
track of this, we also record the step length Ns of bridges, 
denoting the number of linkers n that connect one cluster 
to the next. We find, however, that the distribution of 
the step length decays rapidly with suggesting that 
the initial approximation that most bridges are one-step 
is justified - see Fig. 

As Fig. shows, the cluster size increases with increas¬ 
ing attractive LJ strength 5 . Basically, we are increasing 
the importance of energetic over entropic effects. The 
upper bound of 500 on the cluster size is determined by 
our choice of system - as stated we have 500 total binder 
motifs in each simulation (i.e., 1500 binder beads). A 
system that reaches a cluster size of 500 is therefore fully 
collapsed. For very low values of e, the cluster size is 1 
which means no binder is connected to any other binder, 
and the system is swollen. The cluster size measure is 
thus able to distinguish between the extremes of swollen 
random coil phase and collapsed state. For intermediate 
values of 5 , the system settles into the network phase (see 
also Fig. 1^ and exhibits a finite cluster size, below 500 
but significantly greater than one. 

The other way in which we can alter the balance be¬ 
tween energy and entropy is by increasing the linker 
length. The longer the highly flexible linker chains be¬ 
come, the more configurational entropy they possess. 
This is indeed borne out by Fig. which shows that 
for short linkers, we generally see full collapse but that 
for larger values of Ni^ the system swells and - at a rate 
that depends on e - reverts to the fully random swollen 
state. 

Clearly, the combination of 5 and Ni determines the 
size of the clusters, as well as the connectivity between 
them. We may now ask what the projected mechanical 
properties of the resultant network phase will be. These 
too will depend on the structure of the network: in a di¬ 
lute system, such as the ones we consider here, both the 
random coil state and the collapsed state are expected 
to have poor mechanical performance, if not complete 
lack of rigidity - below the overlap concentration, distinct 
molecules will not entangle to any significant degree in 
the random coil phase, and most certainly will not do so 
in the collapsed state. In the intermediate, network phase 
different molecules will perticipate in a single, connected 
network of flowerlike micelles. In the following, we esti¬ 
mate the mechanical modulus of such a network based 
on the structural features we have just discussed. 



FIG. 9: Average cluster size vs e. Networks with large £ tend 
to form bigger clusters. This tendency is more pronounced in 
molecules with shorter linkers. 



FIG. 10: Average cluster size vs linker length. Increasing 
the proportion of linker beads increases the contribution of 
the entropy of the linker to the overall free energy, favoring 
smaller clusters. 

VI. THE MECHANICAL RESPONSE 

In the simplest approximation, we shall consider the 
effective network that emerges in the regime of interme¬ 
diate 5 and Ni as a rubber with Guassian linkers con¬ 
necting crosslinkers of varying functionality. In previous 
work, Yeo et al. [22] established that the shear modulus 
G of such a network may be computed according to a 
modified classical rubber model (see, for instance, [23|) , 
as 

G = g{f)veknT (3) 

In this equation, z/e is the number concentration of elas¬ 
tically effective chains, and g is the so-called front factor, 
which is to be determined according to the functionality 
distribution. In principle, Uq counts both the contribu¬ 
tions from physical entanglements {vp) and crosslinkers 
(z^c), but since our network - outside of the network phase 
and at sufficiently low concentrations - does not possess 
entanglements that contribute to the modulus we will 
set Ve = ^c- The front factor g{f) accounts for changes 
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in density due to the contraction of the network upon 
crosslinking, as well as for the role of functionality. The 
combined result, obtained first by Yeo and based upon 
earlier work by Duiser and St aver man m, Graessly [25] 
and Tobolsky then yields, for a network of average 
functionality / 


G = 


7-2 

7 


{r^)oJ 


Uc^bT 


(4) 


in this equation is the ratio of the mean 

squared end-to-end distance of the polymer chains in the 
crosslinked network to that of the same chains in the 
uncrosslinked state. Becasue, as we have seen, primarily 
unit step length bridges occur we will count only those as 
elastically effective in the following. Using Eq. we may 
now estimate the shear modulus of the various network 
states in our simulations, using the functionality distri¬ 
bution to compute /, and counting the number density 
of clusters to determine z^c- Sample results are presented 
in Fig. The general trends we observe are consitent 
with what the phase diagram also suggests: As a function 
of both linker length and 5 , there exists and intermediate 
regime where cluster sizes are sufficiently large to permit 
high functionalities (many opportunities for connecting 
to other clusters), but are not yet so large that the num¬ 
ber of potential partner clusters becomes limiting and 
the clusters become single collapsed entities. For shorter 
linker lengths (Fig. [IT] ) the dropoff in modulus at higher 
values of 5 is very pronounced, as the clusters quickly 
become fully isolated - for larger linker lengths (Fig. 0 , 
the wider reach of every cluster (even when it is already 
fairly large) allows the system to retain some connectiv¬ 
ity even at larger cluster sizes. Similar figures may be 
drawn for the dependence on the linker length, and in 
Fig. [^we collect these into a modulus diagram. 

This diagram illustrates what we feel is a crucial point 
about these networks: more is not always better for in¬ 
creased rigidity. An optimum in modulus exists at in¬ 
termediate £ as well as Y/, where the balance between 
connectivity opportunities and functionality is optimal 
for overall mecahnical response. It would be interesting 
to explore whether indeed such an optimal regime exists. 

Finally we note that, obviously, fancier models are 
available to predict the mechanical properties of net¬ 
works such as those we consider. Once fully formed, the 
structure of the hydrogel is no different from that of a 
telechelic gel, which was shown to be well-described by 
the classical theories of Flory m and Stockmayer [28] . 
and whose kinetics of aggregation [29] and de-aggregation 
may be used to assess even their visco-elastic response 

[30l|3l]. 


VII. COMPUTATIONAL RHEOLOGY 

In order to verify the accuracy of rubber elastic mod¬ 
els, we now turn to a direct, quantitative characteriza¬ 
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FIG. 11: Rubber modulus (according to Eq. vs £ for 
B3L12. At large £, the high degree of association of the clus¬ 
ter suppresses the potential for connecting to other clusters - 
large clusters spaced further apart. Thus, the modulus of the 
network drops off at higher e. The solid line is a guide to the 
eye. 



FIG. 12: Overview of computed moduli for all molecular ar¬ 
chitectures studied here. Mechanical performance that is opti¬ 
mal in the sense that it provides the highest modulus arising 
due to a favorable balance between cluster size and cluster 
functionality is obtained for B3L75 and £ = 10. 


tion of the viscoelastic response of the self-assembled 
network. Experimentally, establishing these properties 
is often challenging. There are various rheological meth¬ 
ods, including micro- and macroscopic probing methods 
that can characterize the mechanical behavior of poly¬ 
mer networks. The macroscopic method is called bulk 
rheology. By this method one can relate the shear de¬ 
formation (i.e., the strain 7 ), to the shear response (the 
stress a) in a sample. We copy this protocol and imple¬ 
ment it directly in our simulations. Fig. 13 shows a 














simulation box under oscillatory shear. For purely elas¬ 
tic response this relation is a = Gy, where G is called 
the shear modulus. For viscoelastic materials, that ex¬ 
hibit viscous responses as well as elastic responses, G is 
usually decomposed into a real and an imaginary part: 


FIG. 13: In our oscillatory shear simulations the box is de¬ 
formed (right) in the xy-plane, and executes a periodic oscil¬ 
lation characterized by a frequency uj and an amplitude 70 . 
The network contained inside will accomodate this dynamic 
deformation differently depending on its ability - or inability 
- to undergo structural relaxations on the timescale of the 
applied strain. 


G{uj) = G'{uj)^iG"{uj). (5) 

In this equation, G' is the storage modulus, that char¬ 
acterizes the elastic response of the material. The loss 
modulus G"{uj) quantifies the viscous response. For a 
purely Hookean solid, G'{uj) is constant, and G"{uj) = 
0. For a Newtonian fluid, conversely, G' {uj) = 0 and 
G"{uj) = UJT] with 77 the dynamic viscosity. In oscillatory 
rheology, a time-dependent strain of the form 



0 200 400 600 800 1000 


time 

FIG. 14: The phase lag between sinusoidal signal of deforma¬ 
tion and response determines the elastic and viscous modulus. 
In this figure an example of (5 = J is shown. This is an ex¬ 
ample of an ideal viscoelastic response G' = G". 


cF{t,uj,^o) = 7o (t,cj, 7 o) sin(ncjt) ( 8 ) 

n 

+G^(t, UJ, 7 o ) cos(ncjt} . (9) 

Where G^ is storage modulus, a measure of elastic en¬ 
ergy stored in the material. Provided the applied strain 
is sinusoidal, G^ indicates the magnitude of the har¬ 
monic of the stress response. The loss modulus G^ is 
correlated with the viscous dissipated response and mea¬ 
sures the out-of-phase component of the stress. If one has 
the full signals cr(t) and 7 (t), all the moduli G^ can be 
determined from the complex coefficients of the dis¬ 
crete Fourier transform of the discrete stress time series 
CTt 




7 = 7 osina;t, ( 6 ) 

where 70 , amplitude, is the maximum deformation ap¬ 
plied in each cycle. The time lag between two sinusoidal 
signals determines the viscoelastic moduli of the system: 


G' = — cos S and G" = — sin S. 


7 


7 


(7) 


In the above equations, S is the phase lag between the 
stress and strain signals. Th e phase lag for a viscoelastic 
system is shown in Fig. 


14 


In the case of linear response, this is the entire story 
- stress and strain always have the same frequency as 
no higher harmonics are generated. The only degree of 
freedom that linear response allows is the phase lag S, 
but in linear response the moduli themselves cannot in 
any way be functions of time, or of the strain amplitude, 
themselves. Since we will be interested in the non-linear 
(finite strain) response as well, it is useful to expand the 
stress response in the higher harmonics, expressing it as 
a Fourier series 


1 

= V S (10) 

t=0 

where N is the period of the applied strain times the 
sampling frequency. Using the relationships between the 
Fourier coefficients + C-n and bn = i{cn + c_n), 

along with knowledge that c_^ is the complex conjugate 
of Cn, each harmonic modulus can then be determined as 


/o/ _ 

^ 

7o 


2^(Cn) 

7o 


and G^ 


On 

7o 


2jR(Cn) 

7o 


( 11 ) 


With ^{cn) and 5R(cn) being the imaginary and real 
parts of the complex Fourier coefficients respectively. To 
assess whether or not a particular system, exposed to a 
strain with period uj and amplitude 70 can be consid¬ 
ered in a regime of linear response, we typically monitor 
the ratio of the first to the third harmonic modulus (the 
second is zero, by symmetry). 
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A. Model and simulation protocol 

As stated in section [TTj in all of simulations so far we 
have studied a single, very long chain consisting of 500 
repeats of the binder+linker motif. In the experiments, 
chains are typically much shorter. We expect there to 
be little difference between studying one long chain, or 
chopping this long chain up into many smaller chains. 
Now, we verify this explicitly: we compare the results 
of one repeating chain of 500 binder-linker groups with 
those obtained for 27 individual chains of 19 binder-linker 
groups each (this is the typical repeat number for the ex¬ 
perimentally used Poly(8kU4U) molecule). Both systems 
are allowed to self-assemble into a crosslinked network. 
Each binder is made up of three spherical beads and ev¬ 
ery linker is made up of N beads, and we let N vary 
corresponding to the length of the linker from 3 to 150 
beads. To ensure a fixed bond length in our simulations, 
we choose a large value for the strength of the bond po¬ 
tential ki) = dOO/ceT/cr^. The hydrophobic chains in the 
experimental system are quite short, roughly 5nm, which 
makes them act like a stiff rod [4]. To simulate this we 
choose a large value for the strength of bending potential 
ka = 50 /cbT. We also compare the results of simula¬ 
tions for harmonic bond potential with FENE bond po¬ 
tential [32] by which the nearest-neighbour beads along 
the chain interact through an anharmonic, finitely exten¬ 
sible, non-linear, elastic potential given by 

for r < Rq. 

( 12 ) 

The parameter values of i?o = 1-2 and k = 100 prevent 
entanglement and overlapping of chains [33]. The FENE 
potential diverges logarithmically as r ^ i?o, providing a 
finite distance between chain beads. It is used in coarse 
grained simulations to encode the finite extensibility of 
real polymers. 

According to the above equation, particles closer 
than a interact via a repulsive Lennard-Jones potential 
whereas beyond a the interaction is zero to ensure real 
chain scaling - See also section [TT| We study the effect 
of temperature and density on the mechanical properties 
of the system and compare them with the experimental 
results. To simulate experiments at different concentra¬ 
tions, we keep the number of the beads in simulation 
constant, but change the size of the simulation box to 
decrease or increase the concentration. To do so, we per¬ 
form NET simulations: here, we dial in the pressure to 
control the box volume. We note, that this may induce 
other effects besides those of concentration alone, as we 
necessarily have to go to elevated pressures to realize 
higher concentrations. We have not studied these pre¬ 
stress effects separately. We choose three different initial 
pressures to work with: P = 200, 250 and 300 (molecular 
dynamics units [40]). Our procedure for the equilibration 
and determining the pressure and volume of the system 


UpENEir) = —-kRlln 




is as follows: 



FIG. 15: Left: the system as it is started off, from a random 
initial geometery in a big box. Right: the pressurized system, 
self-assembled in an NPT run 


All simulations start from a random geometry in an 
enormously large box. In the first stage of the simula¬ 
tion, we perform a NVT simulation, that is coupled to 
a Nose -Hoover thermostat, in a large box to thermostat 
the temperature of the system to a desired temperature 
(figure 15). In the next stage, we shrink the simulation 
box to a very high density via a NPT ensemble simu¬ 
lation while also self-assembly of the block copolymers 
takes place. Once the temperature and the total energy 
of the system become stationary, we perform a series of 
NPT runs to expand the box until the pressure is just 
zero. This is the volume at which the chains precisely fit 
inside the box, but do not exert any outward (or inward) 
pressures on it. We take this volume to be the proper 
volume of a particular self-assembled configuration, and 
use it to determine the density. This density, 0*, cor¬ 
responds - though likely not identical - to the overlap 
concentration in experimental systems which is: 


N 

(13) 

rig 

where N is degree of polymerization and Rg is radius 
of gyration of polymer m- Since the number of beads 
in our simulations do not vary, we can say that CjC = 
0/0*. This run is followed by more NPT runs which the 
simulation box is compressed to provide higher density 
systems. 

The rheology simulations are performed with periodic 
boundary conditions using the NVT ensemble, which 
serves to ensure that the pressure of the system changes 
corresponding to the applied deformation - in LAMMPS 
the stress must be determined from the pressure tensor. 
We solve the SLLOD equations of motion which were 
proven to be equivalent to Newton’s equations of motion 
for shear [35]. In this method, instead of boundary driven 
deformation , where shearing deformation is induced to 
the particles by the motion of the boundaries, a velocity 
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gradient is generated to move all the particles propor¬ 
tional to the deformation of the boundaries. These equa¬ 
tions are coupled to a Nose -Hoover chain thermostat in a 
velocity Verlet formulation. The oscillatory shear strain - 
according to Eq. (|^ - is applied to the simulation box in 
the x^-plane, for various amplitudes and oscillation peri¬ 
ods. To obtain the shear stress response, we compute the 
xy component of the symmetric (Cauchy) stress tensor, 
axy, for every bead in the simulation box and sum over 
all atoms every 5 simulation time-steps. In the following 
sections we investigate the feasibility of performing rhe¬ 
ology with the LAMMPS molecular dynamics package 

m- 


B. Computing stress and temperature control 

The oscillatory shear simulation imparts a continu¬ 
ously changing deformation to the simulation box. As 
a result, for affine deformation, each atom (bead) in the 
simulation box can be thought of as being forced to drift 
at a given velocity. For example, if the box is being 
sheared in the x^-plane the atoms at the bottom of the 
box (low Y) have a smaller velocity in the x-direction 
than those atoms at top of the box (at high y). LAMMPS 
subtracts this spurious position-dependent drifting veloc¬ 
ity from each atom while shearing. 

To obtain the viscoelastic modulus of a molecular sys¬ 
tem using Eq. Q, one needs to compute the stress re¬ 
sponse of each particle inside the system to external de¬ 
formation. The Cauchy stress tensor, a 3 by 3 tensor, 
completely defines the state of stress at any particular 
point of a material structure in a given deformed state. 
Written out in components, the Cauchy stress tensor has 
the following form: 


Thus, aij is the force in the ^-direction on the surface 
whose normal is in the j-direction. The diagonal compo¬ 
nents of the stress tensor are force per area in all three 
dimensions, and contain the hydrostatic pressure as well 
as any normal stresses that develop inside the material. 
In our simulations, the complete stress tensor for atom i, 
multiplied by the volume that the atom occupies, in our 
simulations is computed as follows: 


I 

aij = -[mViVj + - '^{ruFij + r2iF2j) 
n=l 

n=l 
^ Na 

+ 2 + 'TZiFzj)]- (15) 



FIG. 17: The stress signal for 7 = 8 %, t = 4 and u: = 0 . 2 , 
reconstructed from only the first harmonic of the full Fourier 
expansion fits very well to the numerical stress data, evidenc¬ 
ing that we am in a regime of linear response. 


CF = 


^xx 
^yx 
^ zx 


^xy 

^yy 

^zy 


^xz 

^yz 

(^ZZ 


(14) 




FIG. 16: Stress fluctuations. Here we show the stress response 
over one period of the strain, for two different maximal shear 
values. The stress response is drawn in red every 50 cycles, 
and the blue line is the average of all the cycles, (a) The shear 
rate is slow compared to fluctuation timescales, so the stress 
response is slightly different in every cycle, (b) For higher 
deformation, larger than 7 8 %, the system becomes more 

rigid and fluctuations are suppressed. 


The first term is a kinetic energy contribution for atom 
i. The second term is a pairwise energy contribution 
where n runs over the Ny neighbors of atom i which are 
all the atoms that are within the cut-off range of the 
potential, t\ and are the position of the two atoms in¬ 
volved in the pairwise interaction, and Fi and F 2 are the 
forces on the 2 atoms resulting from the pairwise inter¬ 
action. The third term is a bond contribution of similar 
form for the W bonds which atom i is part of. There 
is also a similar term for the Na angle interactions that 
atom i is involved in. The so-called virial stress tensor 
is similar to the Cauchy stress tensor, and has all the 
above terms except for the contribution from the kinetic 
energy. In our systems we find no significant difference 
in the moduli calculated by these two different methods 
of computing the stress tensors, because the actual ve¬ 
locities remain low and the energy is dominated by the 
bonded and non-bonded contributions. The com parison 
is presented in Fig. 


20 . As discussed in section VHA 


we convert the stress signal to dynamic moduli using dis¬ 
crete Fourier function. In the linear regime only the first 
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Fourier harmonic of the stress signal has contribution to 
the modulus. To check that the calculated modulus is 
in the linear regime we reconstruct the stress signal from 
only the first Fourier harmonic using the following equa¬ 
tion 


cr(t) = Go ^+7o ^ (G^ sin(27rnt) + G^ cos(27rnt)), 

n=l,3,5,... 

^ ( 16 ) 

where Gq is the n = 1 component of Eq. ( pT] ). The 
reconstructed stress, obtained from the above formula, 
fits very well to the stress data from simulations. This is 
shown in the Fig. 


17 



number of degrees of freedom present in the system: 


-Ekin — 2^BT77/dof 

rj, _ 2-Ekin 

k^Tn^oi 


(17) 


Since the kinetic energy is a function of particle velocity, 
there is often a need to distinguish between a particle’s 
streaming velocity which occurs due to group motion of 
aggregated particles and its thermal velocity due to ther¬ 
mal fluctuation. The sum of the two is the particle’s 
total velocity. Using the Nose -Hoover equations [36] we 
thermostat the translational velocity of particles and sub¬ 
tract a velocity bias that is the result of deforming the 
simulation box. Hence the dependence of the computed 
stress on temperature is only the contribution of thermal 
velocity which is due to the fluctuations of the particles. 



FIG. 18: Strain dependent modulus for B3L45, t = 4, 0/0* = 
2.7. The rise of modulus in the high shear rate region is due to 
bond stretching and hard core potential of the beads during 
the simulation. This regime is unphysical - at such high strain 
rates the implicit solvent interactions break down and our 
simulations can no longer capture the actual deformations 
of the material. We discard points after the minumum, and 
interpret the rapid dropoff as yielding behavior for our system. 


In molecular dynamics simulation of oscillatory shear, 
care should be taken in choosing the right deformation 
rate. If the box deformation rate is larger than the time- 
scale in which the beads interact and fluctuate then the 
position and interaction of the beads are the same in all 
of the oscillatory cycles (see Fig. [^) and the response 
is similar to the response of a crystalline materials which 
has significantly higher modulus than soft materials. Ob¬ 
viously increasing the shear rate eventually would lead 
in stretching the bonds more than equilibrium bond dis¬ 
tance which breaks the bonds and stops the simulation. 
In the Fig. 18 is shown that the mechanical response 
of the system increases for 7 > 8 which is towards sup¬ 
pressing the fluctuations and solid-like structures. This 
rise in modulus eventually terminates at 7 = 25 where 
the simulation stops because of excessive bond stretching. 

In MD simulations a thermostat must be used as a 
means of controlling the temperature of particles. Typi¬ 
cally a target temperature T is specified by the user, and 
the thermostat attempts to equilibrate the system to the 
requested T. To compute the temperature, the kinetic 
energy is divided by the Boltzmann constant and ndof, 


(a) 


(b) 


G' single chain 
-e- G" single chain 

• G' split chain 
-0 - G" split chain 


0.1 0.2 0.5 1.0 2.0 5.0 

Strain [%] 

(c) 

FIG. 19: Simulation snapshots of (a) single long chain with 
500 repeating binder+linker blocks and (b) split chain, where 
the separate chains are distinguished by different colors, (c) 
Strain-dependent simulation results for two different systems 
of single chain (solid line) and split chain (dashed line). Stor¬ 
age modulus is shown in blue color and loss modulus is shown 
in red color. As may be seen here, the results are identical for 
the split and the long system. 



VIII. RESULTS 

Now, we present the results of a series of oscillatory 
rheology simulations of the self-assembling multiblock 
copolymer system. First, we study the effect of changing 
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the molecular architecture from a siugle loug chaiu to sev¬ 
eral short chaius: Straiu-coutrolled simulatious results at 
differeut deusities are showu iu Fig. 19 . To justify the 


validity of our coarse-graiued model to simulate a block 
copolymer, we compare the results of two differeut sys¬ 
tem, (i) a oue repeatiug chaiu that is self-assembled to 
a flower like uetwork {ii) split chaius which are 27 chaius 
of 19 hydrophobic-hydrophilic block each. The system 
sizes iu both cases are equal to 24000 particles. Fig. 
shows equivaleut results for both systems. 

Secoud, we demoustrate the practical equivaleuce (iu 
our simulatiou settiugs) for the virial aud Cauchy stress 
approaches [371 Eg. Here we compare modulus of oue 
system usiug Cauchy aud Virial stresses. As showu iu 
Fig. 


small. 


20 , the differeuce of two stresses iu our model is 



G', Virial 
^ G", Virial 

- G', Cauchy 
-© - G", Cauchy 


0.2 0.5 1.0 2.0 

Strain [%] 


network at 0/0* = 2.7, Cc; = 1. The complex modulus was 
calculated from the two different stress tensors: Cauchy and 
Virial. 



Strain [%] 

(a)Harmonic spring bond, 0/0* = 2.7 



FIG. 21: Strain-dependent comparison between a bead-spring 
polymer network of harmonic bonds vs a bead spring network 
of FENE bonds for B3L45, a; = 1 and T = 2. The storage 
modulus is shown in solid symbols, the loss modulus is shown 
in open symbols. 


Finally, we examine the effects of the details of the 
bond potential. As discussed in section [Vml FENE 
bonds are often used as finitely extensible bonds in coarse 
grained MD simulations to better reflect the finite back¬ 
bone length of polymers. To compare the effect of har¬ 
monic bonds vs. FENE bonds in the shear deformation 
simulations we compare the modulus for two systems con¬ 
taining two above bonds in Fig. 21 . Again, the results 


22 


are completely similar: This may be understood from the 
fact that in the regime where we measure, the non-linear 
stretching regime is never engaged and the FENE bonds 
act as linear springs. 

The effect of concentration is also shown in Fig 
The overlap concentration - which is the concentration 
at which the single molecule network precisely fits within 
the box 0 * - is obtained by relaxing the network volume 
until the pressure first reaches zero. Higher concentra¬ 
tions are shown in units of this overlap concentration. 


The temperature of the simulation also influences the 
modulus of the system. Using Eq. 0 , we can change 
the temperature in the simulation. It affects the mechan¬ 
ics in two ways: first, the self assembly process is affected 
as elevated temperatures favor high entropy states more 
than high energy states, which skews the architectures 
towards smaller cluster sizes. This tends to lower G'. 
Secondly, the mechanical response due to the linkers is 
itself temperature dependent: entropically elastic effects 
scale with temperature. This latter effect would raise 
the storage modulus as a function of temperature 
the case in rubber elasticity. In the Fig. 


23 


as IS 
the strain 

dependent shear modulus for 7 < 20 % at different tem¬ 
peratures is shown. The modulus is calculated based on 
the stress data from simulations that include all the con¬ 
tributions of the stress tensor. The result shows a linear 
regime up to 7 = 5%, after which G' drops off rapidly. 
Linear modulus as a function of temperature and density 
is shown in a logarithmic plot in Fig. 


24 . The slope 


of the power law fit to the temperature dependant lin¬ 
ear modulus from numerical results is 1.28, that is close 
to theoretical prediction for spherical flower-like micelles, 
i.e. 1.95. 
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strain [%] 


Strain [%] 


(a), 0/0* = 2.7 


(b)0/0* = 2.9 



(c)0/0* =3.1 


FIG. 22 : Dynamic moduli G' and G"^ measured in strain-controlled settings at cj = 0.2 and T = 1 at various concentrations. 
The storage modulus is shown in solid symbols, the loss modulus is shown in open symbols. 


23 


IS 


The yield point for all of the systems in Fig. 
where the G' starts to drop, at 7 ~ 5%. That this is 
the yield point is also evidenced by the concurrent and 
characteristic rise in G". To investigate the effect of the 
temperature on the modulus we plot the elastic mod¬ 


ulus for all temperatures as a function of strain in 25 


The experimental measurements at large strains gives a 
quantitative comparison of the contribution of the elas¬ 
tic response at different temperatures [39]. A decrease in 
the magnitude of the modulus with the temperature is 
observed which is in good agreement with the simulation 
results we find here Fig. 25 . The fact that the overall ef¬ 
fect of temperature in these networks is to lower the mod¬ 
ulus suggests that, apparently, the effect of temperature 
on the aggregate size is dominant over the single-chain 
elastic contribution. This suggests that rubber elasticity 
theories should be modified to include also the tempera¬ 
ture dependent changes in connectivity to fully capture 
the mechanical response of the system (and that, there¬ 
fore, the model presented in section |TT| may not offer a 
complete description.) 


where small clusters of the hydrophobic blocks are con¬ 
nected by (typically) single hydrophilic linker chains to 
other such clusters, giving rise to a ’’single molecule net¬ 
works”. The connectivity of this network is generally 
determined by a combination of the size of an aggregate, 
which sets the number of outgoing linkers and thus the 
potential for bridges, and the number of aggregates which 
sets the number of potential partners for such bridg¬ 
ing. At intermediate values of e and Ni, the combina¬ 
tion between ’’supply and demand” of linkers is optimal 
which should result, in the highest moduli for the net¬ 
work material. The multiblock amphiphile system thus 
allows direct control over its supramolecular arrangement 
through molecular design, and thus to mechanical prop¬ 
erties. This suggests much richer applications for these 
systems than what has been established thus far, and in 
particular makes them suitable candidates for exploring 
future biomimetic mechanical performance. 
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